rm(list = ls())

#setwd("YOUR-WORKING-DIRECTORY")

# install.packages(c("doParallel", "foreach", "magic"))

library("doParallel")
library("foreach")
library("magic")
library("MASS")
library("beepr")
library("rr")

set.seed(88)

# # stopCluster(cl)
# cl <- makeCluster(detectCores())
# registerDoParallel(cl)

source("replication-functions.R")


finished <- read.csv("rrreg-fig-1.csv")

# overall.1000 <- subset(finished, case == "overall" & n == 1000)
# overall.2500 <- subset(finished, case == "overall" & n == 2500)
# overall.1e04 <- subset(finished, case == "overall" & n == 1e04)

subgroup.1000 <- subset(finished, case == "subgroup" & n == 1000)
subgroup.2500 <- subset(finished, case == "subgroup" & n == 2500)
subgroup.1e04 <- subset(finished, case == "subgroup" & n == 1e04)

corr1.1000 <- subset(finished, case == "corr1" & n == 1000)
corr1.2500 <- subset(finished, case == "corr1" & n == 2500)
corr1.1e04 <- subset(finished, case == "corr1" & n == 1e04)

corr2.1000 <- subset(finished, case == "corr2" & n == 1000)
corr2.2500 <- subset(finished, case == "corr2" & n == 2500)
corr2.1e04 <- subset(finished, case == "corr2" & n == 1e04)

corr3.1000 <- subset(finished, case == "corr3" & n == 1000)
corr3.2500 <- subset(finished, case == "corr3" & n == 2500)
corr3.1e04 <- subset(finished, case == "corr3" & n == 1e04)


# corrs <- c(1.0, 0.0, 0.0, 
# 	       0.0, 1.0, 0.0, 
# 	       0.0, 0.0, 1.0)


# pop.data <- generateComparisonDataRr(n = 1e07, corrs = corrs, deltas = rep(0.5, 3), p.t = 0.5, n.groups = 5)

# # corrs2 <- c(1.0, 0.0, 0.5, 
# # 	        0.0, 1.0, 0.0, 
# # 	        0.5, 0.0, 1.0)

# # pop.data2 <- generateComparisonDataRr(n = 1e07, corrs = corrs2, deltas = rep(0.5, 3), p.t = 0.5, n.groups = 5)

# # pdf("figure-3.pdf", width = 8.5, height = 5.5)
# # {

# # 	par(mfrow = c(2, 5), mar = c(0, 2, 2, 0), oma = c(2, 11, 2, 2), cex = 0.8)
# # 	x1.h <- hist(pop.data$x1[pop.data$group == 1], plot = FALSE)
# # 	plot(x = x1.h$density, y = x1.h$mids, type = "l", xlim = c(0, 0.5), ylim = c(-2.5, 2.5), 
# # 		xaxt = "n", xlab = "", ylab = "", main = "Group 1")
# # 	abline(h = 0)
# # 	for (i in 2:5) {
# # 		x1.h <- hist(pop.data$x1[pop.data$group == i], plot = FALSE)
# # 		plot(x = x1.h$density, y = x1.h$mids, type = "l", xlim = c(0, 0.5), ylim = c(-2.5, 2.5), 
# # 			xaxt = "n", yaxt = "n", xlab = "", ylab = "", main = paste("Group", i))
# # 		abline(h = 0)
# # 	}

# # 	x1.means.2 <- tapply(pop.data2$x1, pop.data2$group, mean)
# # 	x1.h <- hist(pop.data2$x1[pop.data2$group == 1], plot = FALSE)
# # 	plot(x = x1.h$density, y = x1.h$mids, type = "l", xlim = c(0, 0.5), ylim = c(-2.5, 2.5), 
# # 		xaxt = "n", xlab = "")
# # 	abline(h = 0, lty = 2)
# # 	abline(h = x1.means.2[1])
# # 	for (i in 2:5) {
# # 		x1.h <- hist(pop.data2$x1[pop.data2$group == i], plot = FALSE)
# # 		plot(x = x1.h$density, y = x1.h$mids, type = "l", xlim = c(0, 0.5), ylim = c(-2.5, 2.5), 
# # 			xaxt = "n", yaxt = "n", xlab = "", ylab = "")
# # 		abline(h = 0, lty = 2)
# # 		abline(h = x1.means.2[i])
# # 	}

# # 	mtext("X1", side = 2, outer = TRUE, cex = 0.8, line = 0.5, at = 0.7)
# # 	mtext("X1", side = 2, outer = TRUE, cex = 0.8, line = 0.5, at = 0.2)

# # 	mtext("No Segregation\n on X1", side = 2, outer = TRUE, line = 2.5, at = 0.7, las = 1, font = 2)
# # 	mtext("Segregation\n on X1", side = 2, outer = TRUE, line = 2.5, at = 0.2, las = 1, font = 2)

# # }
# # dev.off()

# necFunc <- c("adiag", "coef.rrreg", "vcov.rrreg", "ginv")

n.bootstrap <- 5000


# # # case 1
# # {

# # 	h <- mean(pop.data$z)

# # 	overall.1000 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# # 	  generateComparisonRr(n.sample = 1000, data = pop.data, h = h, grouping = "overall", weighting = "cue")	

# # 	}


# # 	overall.2500 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# # 	  generateComparisonRr(n.sample = 2500, data = pop.data, 
# # 		h = h, grouping = "overall", weighting = "cue")	

# # 	}


# # 	overall.1e04 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# # 	  generateComparisonRr(n.sample = 1e04, data = pop.data, 
# # 		h = h, grouping = "overall", weighting = "cue")	

# # 	}

# # }


# # case 2
# { 

# 	h <- tapply(pop.data$z, pop.data$group, mean)

# 	subgroup.1000 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 1000, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}


# 	subgroup.2500 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 2500, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}


# 	subgroup.1e04 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 1e04, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# }


# # case 3
# {

# 	corrs <- c(1.0, 0.0, 0.5, 
# 		       0.0, 1.0, 0.0, 
# 		       0.5, 0.0, 1.0)

# 	pop.data <- generateComparisonDataRr(n = 1e07, corrs = corrs, 
# 		deltas = rep(0.5, 3), p.t = 0.5, n.groups = 5)

# 	h <- tapply(pop.data$z, pop.data$group, mean)

# 	corr1.1000 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 1000, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# 	corr1.2500 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 2500, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}


# 	corr1.1e04 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 1e04, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# }


# # case 4
# { 

# 	corrs <- c(1.0, 0.0, 0.5, 
# 		       0.0, 1.0, 0.5, 
# 		       0.5, 0.5, 1.0)

# 	pop.data <- generateComparisonDataRr(n = 1e07, corrs = corrs, 
# 		deltas = rep(0.5, 3), p.t = 0.5, n.groups = 5)

# 	h <- tapply(pop.data$z, pop.data$group, mean)

# 	corr2.1000 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 1000, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# 	corr2.2500 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 2500, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}


# 	corr2.1e04 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 1e04, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# }


# # case 5
# { 

# 	corrs <- c(1.0, 0.5, 0.5, 
# 		       0.5, 1.0, 0.5, 
# 		       0.5, 0.5, 1.0)

# 	pop.data <- generateComparisonDataRr(n = 1e07, corrs = corrs, 
# 		deltas = rep(0.5, 3), p.t = 0.5, n.groups = 5)

# 	h <- tapply(pop.data$z, pop.data$group, mean)

# 	corr3.1000 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 1000, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# 	corr3.2500 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 2500, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}


# 	corr3.1e04 <- foreach(i = 1:n.bootstrap, .combine = "rbind", .export = necFunc) %dopar% {

# 	  generateComparisonRr(n.sample = 1e04, data = pop.data, 
# 		h = h, grouping = "subgroup", weighting = "cue")	

# 	}

# }






pdf(file = "figure-5.pdf", width = 10, height = 10)


{


index <- rep(0:1, n.bootstrap)

par(mfrow = c(3, 3), cex = 0.75, oma = c(2, 2, 2, 12), mar = c(5, 3, 2, 2) + 0.05, xpd = NA)
	# bias
	for (i in 1:3) {

		plot(x = 1:3, y = c(mean(subgroup.1000[index == 0, i]), 
			mean(subgroup.2500[index == 0, i]), 
			mean(subgroup.1e04[index == 0, i])
			), 
	        lwd = 1.5, 
		    type = "b", 
		    ylim = c(0.485, 0.525), 
		    yaxt = "n",
		    ylab = "Bias", 
		    xaxt = "n", 
		    xlab = "Sample Size", 
		    main = substitute(delta[m], list(m = i - 1)), 
		    lty = 1, 
		    pch = 1
		)

		par(xpd = FALSE); abline(h = 0.5, col = "darkgray"); par(xpd = NA)
	    axis(side = 2, at = seq(0.49, 0.52, 0.01), labels = seq(0.49, 0.52, 0.01) - .5)
	    axis(side = 1, at = 1:3, labels = c("1000", "2500", "10000"))

		lines(x = 1:3, y = c(mean(subgroup.1000[index == 1, i]), 
			mean(subgroup.2500[index == 1, i]), 
			mean(subgroup.1e04[index == 1, i])
			), 
		    lwd = 1.5, 
		    type = "b", 
		    col = "black", 
		    lty = 2, 
		    pch = 16
		   )

		lines(x = 1:3, y = c(mean(corr1.1000[index == 1, i]), 
			mean(corr1.2500[index == 1, i]), 
			mean(corr1.1e04[index == 1, i])
			), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "red", 
		    lty = 1, 
		    pch = 2
		   )

		lines(x = 1:3, y = c(mean(corr2.1000[index == 1, i]), 
			mean(corr2.2500[index == 1, i]), 
			mean(corr2.1e04[index == 1, i])
			), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "blue", 
		    lty = 2, 
		    pch = 17
		   )

		lines(x = 1:3, y = c(mean(corr3.1000[index == 1, i]), 
			mean(corr3.2500[index == 1, i]), 
			mean(corr3.1e04[index == 1, i])
			), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "darkgreen", 
		    lty = 1, 
		    pch = 4
		   )

	}


	legend("topright", legend = c("No aux. info",  
	  "No segregation", "Segregation on X1", "Segregation on X1, X2", "X1, X2 correlated"), 
	  col = c("black", "black", "red", "blue", "darkgreen"), 
	  lty = c(1, 2, 1, 2, 1),
	  pch = c(1, 16, 2, 17, 4), 
	  lwd = 1.5, 
	  bty = "n", 
	  seg.len = 1.5, 
	  inset = c(-1, 0)
	  )


	# rmse
	for (i in 1:3) {

		plot(x = 1:3, y = c(rmse(subgroup.1000[index == 0, i], 0.5), 
			rmse(subgroup.2500[index == 0, i], 0.5), 
			rmse(subgroup.1e04[index == 0, i], 0.5)
			), 
		    lwd = 1.5, 
		    type = "b", 
		    ylim = c(0, 0.2), 
		    ylab = "RMSE", 
		    xaxt = "n", 
		    xlab = "Sample Size", 
		    main = substitute(delta[m], list(m = i - 1)), 
		    lty = 1, 
		    pch = 1
		)
	    axis(side = 1, at = 1:3, labels = c("1000", "2500", "10000"))

		par(xpd = FALSE); abline(h = 0, col = "darkgray"); par(xpd = NA)

		lines(x = 1:3, y = c(rmse(subgroup.1000[index == 1, i], 0.5), 
			rmse(subgroup.2500[index == 1, i], 0.5), 
			rmse(subgroup.1e04[index == 1, i], 0.5)
			), 
		    type = "b", 
	        lwd = 1.5, 
		    col = "black", 
		    lty = 2, 
		    pch = 16
		   )

		lines(x = 1:3, y = c(rmse(corr1.1000[index == 1, i], 0.5), 
			rmse(corr1.2500[index == 1, i], 0.5), 
			rmse(corr1.1e04[index == 1, i], 0.5)
			), 
		    type = "b", 
	        lwd = 1.5, 
		    col = "red", 
		    lty = 1, 
		    pch = 2
		   )

		lines(x = 1:3, y = c(rmse(corr2.1000[index == 1, i], 0.5), 
			rmse(corr2.2500[index == 1, i], 0.5), 
			rmse(corr2.1e04[index == 1, i], 0.5)
			), 
		    type = "b", 
	        lwd = 1.5, 
		    col = "blue", 
		    lty = 2, 
		    pch = 17
		   )

		lines(x = 1:3, y = c(rmse(corr3.1000[index == 1, i], 0.5), 
			rmse(corr3.2500[index == 1, i], 0.5), 
			rmse(corr3.1e04[index == 1, i], 0.5)
			), 
		    type = "b", 
	        lwd = 1.5, 
		    col = "darkgreen", 
		    lty = 1, 
		    pch = 4
		   )

	}

	legend("topright", legend = c("No aux. info",  
	  "No segregation", "Segregation on X1", "Segregation on X1, X2", "X1, X2 correlated"), 
	  col = c("black", "black", "red", "blue", "darkgreen"), 
	  lty = c(1, 2, 1, 2, 1),
	  pch = c(1, 16, 2, 17, 4), 
	  lwd = 1.5, 
	  bty = "n", 
	  seg.len = 1.5, 
	  inset = c(-1, 0)
	  )

	# coverage
	for (i in 1:3) {

		plot(x = 1:3, y = c(
				cover(x   = subgroup.1000[index == 0, i], 
					  se  = subgroup.1000[index == 0, i + 3], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = subgroup.2500[index == 0, i], 
					  se  = subgroup.2500[index == 0, i + 3], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = subgroup.1e04[index == 0, i], 
					  se  = subgroup.1e04[index == 0, i + 3], 
					  par = 0.5, 
					  alpha = 0.05)
				), 
		    lwd = 1.5, 
		    type = "b", 
		    ylim = c(0.9, 1), 
		    ylab = "Coverage", 
		    xaxt = "n", 
		    xlab = "Sample Size", 
		    main = substitute(delta[m], list(m = i - 1)), 
		    lty = 1, 
		    pch = 1
		)
		
	    axis(side = 1, at = 1:3, labels = c("1000", "2500", "10000"))
		par(xpd = FALSE); abline(h = 0.95, col = "darkgray"); par(xpd = NA)

		lines(x = 1:3, y = c(
				cover(x   = subgroup.1000[index == 1, i], 
					  se  = subgroup.1000[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = subgroup.2500[index == 1, i], 
					  se  = subgroup.2500[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = subgroup.1e04[index == 1, i], 
					  se  = subgroup.1e04[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05)
				), 
	        lwd = 1.5, 
		    type = "b", 
		    lty = 2, 
		    pch = 16
		   )

		lines(x = 1:3, y = c(
				cover(x   = corr1.1000[index == 1, i], 
					  se  = corr1.1000[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr1.2500[index == 1, i], 
					  se  = corr1.2500[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr1.1e04[index == 1, i], 
					  se  = corr1.1e04[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05)
				), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "red", 
		    lty = 1, 
		    pch = 2
		   )

		lines(x = 1:3, y = c(
				cover(x   = corr2.1000[index == 1, i], 
					  se  = corr2.1000[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr2.2500[index == 1, i], 
					  se  = corr2.2500[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr2.1e04[index == 1, i], 
					  se  = corr2.1e04[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05)
				), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "blue", 
		    lty = 2, 
		    pch = 17
		   )


		lines(x = 1:3, y = c(
				cover(x   = corr3.1000[index == 1, i], 
					  se  = corr3.1000[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr3.2500[index == 1, i], 
					  se  = corr3.2500[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05), 
				cover(x   = corr3.1e04[index == 1, i], 
					  se  = corr3.1e04[index == 1, i + 3], 
					  par = 0.5, 
					  alpha = 0.05)
				), 
	        lwd = 1.5, 
		    type = "b", 
		    col = "darkgreen", 
		    lty = 1, 
		    pch = 4
		   )

		}

		legend("topright", legend = c("No aux. info",  
		  "No segregation", "Segregation on X1", "Segregation on X1, X2", "X1, X2 correlated"), 
		  col = c("black", "black", "red", "blue", "darkgreen"), 
		  lty = c(1, 2, 1, 2, 1),
		  pch = c(1, 16, 2, 17, 4), 
		  lwd = 1.5, 
		  bty = "n", 
		  seg.len = 1.5, 
		  inset = c(-1, 0)
		  )

  # title(main = c("Statistical Properties of the Augmented List Experiment Estimator"), outer = TRUE)


    

}

dev.off()




# finished <- rbind(
# 	# cbind(overall.1000, n = 1000, case = "overall"), 
# 	# cbind(overall.2500, n = 2500, case = "overall"), 
#  #    cbind(overall.1e04, n = 1e04, case = "overall"), 

# 	cbind(subgroup.1000, n = 1000, case = "subgroup"), 
# 	cbind(subgroup.2500, n = 2500, case = "subgroup"), 
# 	cbind(subgroup.1e04, n = 1e04, case = "subgroup"), 

# 	cbind(corr1.1000, n = 1000, case = "corr1"), 
# 	cbind(corr1.2500, n = 2500, case = "corr1"), 
# 	cbind(corr1.1e04, n = 1e04, case = "corr1"), 

# 	cbind(corr2.1000, n = 1000, case = "corr2"), 
# 	cbind(corr2.2500, n = 2500, case = "corr2"), 
# 	cbind(corr2.1e04, n = 1e04, case = "corr2"), 

# 	cbind(corr3.1000, n = 1000, case = "corr3"), 
# 	cbind(corr3.2500, n = 2500, case = "corr3"), 
# 	cbind(corr3.1e04, n = 1e04, case = "corr3")
# 	)

# finished <- as.data.frame(finished)
# names(finished)[7] <- "augmented"

# write.csv(finished, "rrreg-fig-1.csv", row.names = FALSE)


# stopCluster(cl)




# q("no")